example(pcountOpen)

re <- ranef(m1)

plot(re, layout=c(5,5), subset = site %in% 1:25 & year %in% 1:2,
     xlim=c(-1,15))
     
G           # Matrix of recruits from simulation in example


sim_conditional_recruits <- function(x, fm) {
  N <- x ## To be clear that x is abundance
  dynamics <- fm@dynamics
  n_sites <- nrow(x)
  n_years <- ncol(x)
  gam_vec <- predict(fm, type="gamma")$Predicted
  gam <- matrix(gam_vec, n_sites, n_years-1)
  omega_vec <- predict(fm, type="omega")$Predicted
  omega <- matrix(omega_vec, n_sites, n_years-1)
  S <- matrix(NA_integer_, n_sites, n_years-1) # Survivors
  G <- matrix(NA_integer_, n_sites, n_years-1) # Recruits
  rcat <- function(probs) which(rmultinom(n=1, size=1, prob=probs)==1)
  for(i in 1:n_sites) {
    for(t in 2:n_years) {
      possible_survivors <- 0:N[i,t-1]
      possible_recruits <- N[i,t]-possible_survivors
      p_survivors <- dbinom(possible_survivors, N[i,t-1], omega[i,t-1])
      if(dynamics=="constant")
        p_recruits <- dpois(possible_recruits, gam[i,t-1])
      else if(dynamics=="autoreg")
        p_recruits <- dpois(possible_recruits, gam[i,t-1]*N[i,t-1])
      else
        stop("dynamics not implemented")
      p_joint0 <- p_survivors*p_recruits
      p_joint <- p_joint0/sum(p_joint0)
      index <- rcat(p_joint)
      S[i,t-1] <- possible_survivors[index]
      G[i,t-1] <- possible_recruits[index]
    }
  }
  return(G=colSums(G))
}


G_post <- predict(re, func=sim_conditional_recruits, nsims=10, fm=m1)

rowMeans(G_post) ## Estimates of *realized* recruits per time interval
apply(G_post, 1, quantile, probs=c(0.025, 0.975))


colSums(G)  # Actual recruits per time interval (from data simulation)
gam*M       # Estimate of the expected number of recruits per interval